#include <stdio.h>
#include <stdlib.h>
#include "iam_matop.h"

void printmat_F64(IAM_RasterView_F64 mat)
{
    for (int y = 0; y < mat->m_size_y; ++y) {
        for (int x = 0; x < mat->m_size_x; ++x) {
            double val = IAM_RasterView_getValue_F64(mat, x, y);
            printf("%.8f, ", val);
        }
        printf("\n");
    }
}

void testSVD_F64(double* pa, size_t Nx, size_t Ny)
{
    IAM_RasterView_F64 A = IAM_RasterView_init_F64(Nx, Ny, pa);

    printf("A = \n");
    printmat_F64(A);

    IAM_RasterData_F64 U = IAM_RasterData_init0_F64();
    IAM_RasterData_F64 S = IAM_RasterData_init0_F64();
    IAM_RasterData_F64 V = IAM_RasterData_init0_F64();

    IAM_MatOp_svd_F64(A, U, S, V);

    printf("U = \n");
    printmat_F64(U->m_view);

    printf("S = \n");
    printmat_F64(S->m_view);

    printf("V = \n");
    printmat_F64(V->m_view);

    IAM_RasterView_release_F64(&A);
    IAM_RasterData_release_F64(&U);
    IAM_RasterData_release_F64(&S);
    IAM_RasterData_release_F64(&V);
    printf("\n");
}


int main(int argc, char** argv)
{
    double p1[] = {1, 2, 4, 3, 5, 5, 7, 1};
    printf("======== Test 1 ========\n");
    testSVD_F64(p1, 2, 4);
    
    double p2[] = {1, 2, 4, 3, 5, 5, 7, 1, 2};
    printf("======== Test 2 ========\n");
    testSVD_F64(p2, 3, 3);
    
    double p3[] = {1, 2, 4, 3, 5, 5, 7, 1};
    printf("======== Test 3 ========\n");
    testSVD_F64(p3, 4, 2);
    
    return 0;
}

